
#include "visual.h"


void output_uvp(const char *szProblem,
		 double xlength,
                 double ylength,
                 int  imax,
                 int  jmax,
		 double dx,
		 double dy,
                 double **U,
                 double **V,
                 double **P,
		 int **Flag)
{
    static int nCount = 0;  /* counts the calls of this function */
    int i, j;
    double **tmp = matrix( 1, imax, 1, jmax );
    double **Psi;
    double **Zeta;
    double **fluid;
    char szFileName[80];

    printf("output is written.\n");
    if ( nCount==0 ){               /* coordinates for AVS: only once */
      for( i = 1; i <= imax; i++){
        for( j = 1; j <= jmax; j++){
          tmp[i][j] = ((float)i-0.5)*dx;
        }
      }
      sprintf( szFileName, "%s.X.koo", szProblem );
      write_matrix( szFileName,tmp,1,imax,1,jmax,xlength,ylength,nCount==0);
      for( i = 1; i <= imax; i++){
        for( j = 1; j <= jmax; j++){
          tmp[i][j] = ((float)j-0.5)*dy;
        }
      }
      sprintf( szFileName, "%s.Y.koo", szProblem );
      write_matrix( szFileName,tmp,1,imax,1,jmax,xlength,ylength,nCount==0);
    }

    /* U-Matrix */
    for( i = 1; i <= imax; i++){
      for( j = 1; j <= jmax; j++){
        if ( Flag[i][j]&C_F ) tmp[i][j] = (U[i-1][j] + U[i][j]) * 0.5;
        else tmp[i][j] = 0.;
      }
    }

    sprintf( szFileName, "%s.U.out", szProblem );

    write_matrix( szFileName,tmp,1,imax,1,jmax,xlength,ylength,nCount==0);

    /* V-Matrix */
    for( i = 1; i <= imax; i++){
      for( j = 1; j <= jmax; j++){
        if ( Flag[i][j]&C_F ) tmp[i][j]  = (V[i][j-1] + V[i][j]) * 0.5;
        else tmp[i][j] = 0.;
      }
    }

    sprintf( szFileName, "%s.V.out", szProblem );

    write_matrix( szFileName,tmp,1,imax,1,jmax,xlength,ylength,nCount==0);

    /* P-Matrix */
    sprintf( szFileName, "%s.P.out", szProblem );

    write_matrix( szFileName,P,1,imax,1,jmax,xlength,ylength,nCount==0);

    free_matrix( tmp,1,imax,1,jmax );

    /* Fluid Matrix*/
    fluid = matrix(1, imax, 1, jmax);
    init_matrix( fluid, 1, imax, 1, jmax, 0.0);

    for ( i = 1; i < (imax + 1); i++ ) {
	for ( j = 1; j < (jmax + 1); j++ ) {
	    if ( is_fluid(Flag[i][j]) ) {
		fluid[i][j] = 1.0;
	    }
	}
    }

    sprintf( szFileName, "%s.fluid.out", szProblem );

    write_matrix( szFileName, fluid, 1, imax, 1, jmax, xlength, ylength, nCount == 0 );

    free_matrix( fluid, 1, imax, 1, jmax );

    /* write field file to be read by AVS */
    write_fieldfile( szProblem,imax,jmax,xlength,ylength,nCount+1);

    Psi = matrix(0,imax,0,jmax);
    Zeta = matrix(1,imax,1,jmax);

    calculate_psi_zeta(dx,dy,imax,jmax,U,V,Psi,Zeta,Flag);

    sprintf( szFileName, "%s.psi.out", szProblem );
    write_matrix( szFileName,Psi,1,imax,1,jmax,xlength,ylength,nCount==0);

    sprintf( szFileName, "%s.zeta.out", szProblem );
    write_matrix( szFileName,Zeta,1,imax,1,jmax,xlength,ylength,nCount==0);

    free_matrix( Zeta, 1, imax, 1, jmax );
    free_matrix( Psi, 0, imax, 0, jmax );

    ++nCount;
}


void write_fieldfile( 	const char *szProblem,
			int imax, int jmax,
			double xlength,
			double ylength,
			int nCount
			)
{
	char szFileName[80];
	FILE *fp=NULL;

	sprintf( szFileName, "%s.fld", szProblem );

	fp = fopen( szFileName, "w");
	fprintf(fp,"# AVS field file %s\n",szProblem);
	fprintf(fp,"#\n");
	fprintf(fp,"nstep   = %d\n",nCount);
	fprintf(fp,"ndim    = 2\n");
	fprintf(fp,"dim1    = %d\n",imax);
	fprintf(fp,"dim2    = %d\n",jmax);
	fprintf(fp,"nspace  = 2\n");
	fprintf(fp,"veclen  = 6\n");
	fprintf(fp,"data    = float\n");
	fprintf(fp,"field   = irregular\n");
	fprintf(fp,"min_ext = %f %f\n",0.0,0.0);
	fprintf(fp,"max_ext = %f %f\n",xlength,ylength);
	sprintf( szFileName, "%s.X.koo", szProblem );
	fprintf(fp,"coord 1 file=%s filetype=binary skip=0 offset=0 stride=1 close=1\n",szFileName);
	sprintf( szFileName, "%s.Y.koo", szProblem );
	fprintf(fp,"coord 2 file=%s filetype=binary skip=0 offset=1 stride=1 close=1\n",szFileName);
	fprintf(fp,"DO\n");
	sprintf( szFileName, "%s.U.out", szProblem );
	fprintf(fp,"variable 1 file=%s filetype=binary skip=0 offset=0 stride=1 close=0\n",szFileName);
	sprintf( szFileName, "%s.V.out", szProblem );
	fprintf(fp,"variable 2 file=%s filetype=binary skip=0 offset=0 stride=1 close=0\n",szFileName);
	sprintf( szFileName, "%s.P.out", szProblem );
	fprintf(fp,"variable 3 file=%s filetype=binary skip=0 offset=0 stride=1 close=0\n",szFileName);
	sprintf( szFileName, "%s.psi.out", szProblem );
	fprintf(fp,"variable 4 file=%s filetype=binary skip=0 offset=0 stride=1 close=0\n",szFileName);
	sprintf( szFileName, "%s.zeta.out", szProblem );
	fprintf(fp,"variable 5 file=%s filetype=binary skip=0 offset=0 stride=1 close=0\n",szFileName);
	sprintf( szFileName, "%s.fluid.out", szProblem );
	fprintf(fp,"variable 6 file=%s filetype=binary skip=0 offset=0 stride=1 close=0\n",szFileName);
	fprintf(fp,"EOT\n");
	fprintf(fp,"ENDDO\n");
	fclose(fp);
}


void write_particlefieldfile(char *outputfile, int partmax, int nCount){

  FILE *fp=NULL;
  char szFileName[80];

  sprintf( szFileName, "%s.fld", outputfile);

  fp = fopen( szFileName,"w");
  fprintf(fp,"# AVS Partikel-Fieldfile %s\n",outputfile);
  fprintf(fp,"#\n");
  fprintf(fp,"nstep=%d\n",nCount);
  fprintf(fp,"ndim=2\n");
  fprintf(fp,"dim1=%d\n",partmax-1);
  fprintf(fp,"dim2=1\n");
  fprintf(fp,"nspace=2\n");
  fprintf(fp,"veclen=0\n");
  fprintf(fp,"data=float\n");
  fprintf(fp,"field=irregular\n");
  fprintf(fp,"DO\n");
  sprintf( szFileName, "%s.prt", outputfile);
  fprintf(fp,"coord 1 file=%s filetype=ascii skip=0 offset=0 stride=2 close=0\n",szFileName);
  fprintf(fp,"coord 2 file=%s filetype=ascii skip=0 offset=1 stride=2 close=0\n",szFileName);
  fprintf(fp,"EOT\n");
  fprintf(fp,"ENDDO\n");

  fclose(fp);

  ++nCount;
}


struct particle_line *set_particle(	int N,
					double x1,
					double y1,
					double x2,
					double y2
					)
{
	struct particle_line *p = NULL;

	double dx = 0.0;
	double dy = 0.0;

	double x_position = x1;
	double y_position = y1;

	int length = 1;
	int i;

	/* Make sure that we have some values that make sense! */
	/*
	assert ( x2 > x1 );
	assert ( y2 > y1 );
	*/

	dx = (x2 - x1) / N;
	dy = (y2 - y1) / N;

	/* Allocate the space for the whole list */
	p = (struct particle_line *) malloc( sizeof(struct particle_line) * N );

	for ( i = 0; i < N; i++) {
		p[i].Particle = (struct particle *) malloc( sizeof( struct particle) );
		p[i].length = length;
		p[i].Particle -> x = x_position;
		p[i].Particle -> y = y_position;

		/* Allocate space for the next item */
		p[i].Particle -> next = (struct particle *) malloc (sizeof(struct particle) );
		p[i].Particle -> next -> x = x_position;
		p[i].Particle -> next -> y = y_position;
		p[i].Particle -> next -> next = NULL;

		x_position += dx;
		y_position += dy;

		/*
		printf("x = %f \t y = %f\n", p[i].Particle -> x, p[i].Particle -> y );
		printf("x1 = %f \t y1 = %f\n", p[i].Particle -> next -> x, p[i].Particle -> next -> y );
		*/
	}
	return p;
}


void move_particle(	struct particle_line *Particles_List,
			int N,
			double dt,
			double dx,
			double dy,
			int imax,
			int jmax,
			double **U,
			double **V,
			int **Flag
		)
{
	int itemp, jtemp, k;
	double xtemp;
	double ytemp;
	double x_1, x_2;
	double y_1, y_2;
	double u_1, u_2, u_3, u_4;
	double u_val;
	double v_1, v_2, v_3, v_4;
	double v_val;
	struct particle *current;
	struct particle *previous;
	struct particle *next;

	itemp = 0;
	jtemp = 0;

	for ( k = 0; k < N; k++ ) {
		previous = NULL;
		current = Particles_List[k].Particle;
		if(current!=NULL) {
				next = Particles_List[k].Particle->next;
		}

		while(current!=NULL) {
			xtemp=current->x;
			ytemp=current->y;

			itemp=(int) (xtemp/dx) + 1;
			jtemp=(int) ((ytemp + dy/2)/dy) + 1;

			x_1=(itemp-1.0)*dx;
			x_2=itemp*dx;
			y_1=((jtemp-1.0)-0.5)*dy;
			y_2=(jtemp-0.5)*dy;

			u_1=U[itemp-1][jtemp-1];
			u_2=U[itemp][jtemp-1];
			u_3=U[itemp-1][jtemp];
			u_4=U[itemp][jtemp];

			u_val=(1/(dx*dy))*((x_2-xtemp)*(y_2-ytemp)*u_1+(xtemp-x_1)*(y_2-ytemp)*u_2+
								 (x_2-xtemp)*(ytemp-y_1)*u_3+(xtemp-x_1)*(ytemp-y_1)*u_4);


			itemp=(int) ((xtemp + dx/2)/dx) + 1;
			jtemp=(int) (ytemp/dy) + 1;

			x_1=((itemp-1)-0.5)*dx;
			x_2=(itemp-0.5)*dx;
			y_1=(jtemp-1)*dy;
			y_2=jtemp*dy;

			v_1=V[itemp-1][jtemp-1];
			v_2=V[itemp][jtemp-1];
			v_3=V[itemp-1][jtemp];
			v_4=V[itemp][jtemp];

			v_val=(1/(dx*dy))*((x_2-xtemp)*(y_2-ytemp)*v_1+(xtemp-x_1)*(y_2-ytemp)*v_2+
							 (x_2-xtemp)*(ytemp-y_1)*v_3+(xtemp-x_1)*(ytemp-y_1)*v_4);

			xtemp = xtemp + dt*u_val;
			ytemp = ytemp + dt*v_val;


			if(xtemp<0.0) {
				itemp=(int) (xtemp/dx);
			}
			else {
				itemp=(int) (xtemp/dx)+1;
			}

			if(ytemp<0.0) {
				jtemp=(int) (ytemp/dy);
			}
			else {
				jtemp=(int) (ytemp/dy)+1;
			}

			if(is_obstacle(Flag[itemp][jtemp])) {
				if(Particles_List[k].Particle==current) {
					Particles_List[k].Particle = next;
					Particles_List[k].length--;
					free(current);
					current = next;
				}
				else {
					previous->next = next;
					Particles_List[k].length--;
					free(current);
					current = next;
				}
			}
			else {
				current->x=xtemp;
				current->y=ytemp;

				previous = current;
				current = next;
			}

			if(current!=NULL ) {
				next = next->next;
			}
		}
	}
}



void insert_particles(	struct particle_line *Particle_Lines, int N)
{
	int i;
	struct particle *p = NULL;
	struct particle *p_new = NULL;

	for ( i = 0; i < N; i++ ) {
		p = Particle_Lines[i].Particle;
		p_new = (struct particle *) malloc (sizeof(struct particle));
		p_new -> x = p -> x;
		p_new -> y = p -> y;
		p_new -> next = p;
		Particle_Lines[i].Particle = p_new;
		Particle_Lines[i].length ++;
	}
}

void write_particles(	struct particle_line *Particle_Lines,
			int N,
			int partmax,
			char *outputfile
			)
{
	FILE * fp = NULL;
	int i = 0;
	struct particle *p = NULL;
	int number_of_particles = 0;
	char prt_file[32];

	sprintf ( prt_file, "%s.prt", outputfile);
	fp = fopen(prt_file, "a");

	if ( fp == NULL ) {
		fprintf(stderr, "Error opening the outputfile: %s\n", prt_file);
		exit(-1);
	}

	/* printf("%d \t %d \n", N, partmax); */
	for ( i = 0; i < N; i++ ) {
		p = Particle_Lines[i].Particle;

		while ( p != NULL /*&& p -> next != NULL*/ ) {

			fprintf(fp, "%f\t%f\n", p -> x, p -> y);
			p = p -> next;
			number_of_particles ++;
		}
		/* TODO: This check is inside the while-loop or the for-loop ?? */

		if ( number_of_particles > partmax ) {
			break;
		}

	}

	/* flush all the data and close the file */
	fflush(fp);
	if ( fclose(fp) ) {
		fprintf(stderr,"Error closing the outputfile: %s\n",prt_file);
		exit( -1);
	}

	write_particlefieldfile(outputfile, partmax, nCount);
	nCount++;

}

void pathlines(	struct particle_line *Particle_Lines,
		int N,
		char *outputfile,
		double dt_path,
		double t,
		double dx,
		double dy,
		double dt,
		int imax,
		int jmax,
		double **U,
		double **V,
		int **Flag
		)
{
	FILE *fp = NULL;
	char prt_file[32];

	sprintf ( prt_file, "%s.prt", outputfile);

	move_particle (
			Particle_Lines,
			N,
			dt,
			dx,
			dy,
			imax,
			jmax,
			U,
			V,
			Flag
		);


	/* create the file in the 1st time step */
	if ( fabs(t - 0.0) < TOL ) {
		fp = fopen(prt_file, "w");
		if ( fp != NULL ) {
			fprintf(fp,"%f\t%f\t%f\n", (float)(imax*dx), (float)(jmax*dy), (float)N);
			fclose (fp);
		}
		else {
			fprintf (stderr, "Error while opening the output file: %s", prt_file);
		}
	}
	else {
		/* do nothing */
	}

	if ( fabs (fmod(t, dt_path) - 0.0) < TOL ) {
		write_particles( Particle_Lines, N, 2*N, outputfile );
	}

}

void streaklines(	struct particle_line *Particle_Lines,
			int N,
			double dt_inject,
			double dt_streak,
			double t,
			double dx,
			double dy,
			double dt,
			int imax,
			int jmax,
			double **U,
			double **V,
			char *outputfile,
			int **Flag
		)
{
	FILE *fp = NULL;
	char prt_file[32];

	sprintf ( prt_file, "%s.prt", outputfile);

	move_particle (
			Particle_Lines,
			N,
			dt,
			dx,
			dy,
			imax,
			jmax,
			U,
			V,
			Flag
		);

	if ( fabs (fmod(t, dt_inject) - 0.0) < TOL ) {
		insert_particles(Particle_Lines, N);
	}

	/* create the file in the 1st time step */
	if ( fabs(t - 0.0) < TOL ) {
		fp = fopen(prt_file, "w");
		if ( fp != NULL ) {
			fprintf(fp,"%f\t%f\t%f\n", (float)(imax*dx), (float)(jmax*dy), (float)N);
			fclose (fp);
		}
		else {
			fprintf (stderr, "Error while opening the output file: %s", prt_file);
		}
	}
	else {
		/* do nothing */
	}

	if ( fabs (fmod(t, dt_streak) - 0.0) < TOL ) {
		write_particles( Particle_Lines, N, STREAK_LINES_PART_MAX, outputfile );
	}

}


void calculate_psi_zeta(	double dx,
				double dy,
				int imax,
				int jmax,
				double **U,
				double **V,
				double **Psi,
				double **Zeta,
				int **Flag
			)
{
	int i, j;
	/* Fill both matrices with zeros initially */
	init_matrix(Psi, 0, imax, 0, jmax, 0.0);
	init_matrix(Zeta, 1, imax, 1, jmax, 0.0);

	for ( i = 1; i <= imax; i++ ) {
		Psi[i][0] = Psi[i - 1][0] - V[i][0] * dx;
	}

	for ( i = 0; i <= imax; i++ ) {
		for ( j = 1; j <= jmax; j++ ) {
			if ( is_fluid(Flag[i][j]) ) {
				Psi[i][j] = Psi[i][j - 1] + U[i][j] * dy;
			}
			else {
				Psi[i][j] = Psi[i][j - 1];
			}
		}
	}

	for ( i = 1; i < imax; i++ ) {
		for ( j = 1; j < jmax; j++) {
			if ( is_fluid (Flag[i][j]) ) {
				Zeta[i][j] = ( U[i][j + 1] - U[i][j] ) / dy - ( V[i + 1][j] - V[i][j]) / dx ;
			}
		}
	}
}


void trace_fluid(
		struct particle_line *Particle_Lines,
		char 	*outputfile,
		int 	N,
		int 	npart,
		double 	t,
		int 	xlength,
		int 	ylength,
		double 	dt_visual
		)
{
	FILE *fp = NULL;
	char prt_file[32];

	sprintf ( prt_file, "%s.prt", outputfile);

	/* create the file in the 1st time step */
	if ( fabs(t - 0.0) < TOL ) {
		fp = fopen(prt_file, "w");
		if ( fp != NULL ) {
			fprintf(fp, "%f\t%f\n", (float)xlength, (float)ylength);
			fclose (fp);
		}
		else {
			fprintf (stderr, "Error while opening the output file: %s", prt_file);
		}
	}
	else {
		/* do nothing */
	}

/*	if ( fabs (fmod(t, dt_visual) - 0.0) < TOL ) { */
		write_particles( Particle_Lines, N, npart, outputfile );
	/*}*/

}
